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In this article we present and compare two different, coarse-grained approaches to model electrostatic in¬ 
teractions of disc-shaped aromatic molecules, specifically coronene. Our study builds on previous work [J. 
Chem. Phys. 141 , 214110 (2014)] where we proposed, based on a systematic coarse-graining procedure 
starting from the atomistic level, an anisotropic effective (Gay-Berne-like) potential capable of describing 
van-der-Waals contributions to the interaction energy. To take into account electrostatics, we introduce, first, 
a linear quadrupole moment along the symmetry axis of the coronene disc. The second approach takes into 
account the fact that the partial charges within the molecules are distributed in a ring-like fashion. We then 
reparametrize the effective Gay-Berne-like potential such that it matches, at short distances, the ring-ring 
potential. To investigate the validity of these two approaches, we perform many-particle Molecular Dynam¬ 
ics (MD) simulations, focusing on the crystalline phase (karpatite) where electrostatic interaction effects are 
expected to be particularly relevant for the formation of tilted stacked columns. Specifically, we investigate 
various structural parameters as well as the melting transition. We find that the second approach yields 
consistent results with those from experiments despite the fact that the underlying potential decays with the 
wrong distance dependence at large molecule separations. Our strategy can be transferred to a broader class 
of molecules, such as benzene or hexabenzocoronene. 


I. INTRODUCTION 

In the present article we analyze two methods to ef¬ 
fectively incorporate electrostatics in symmetric disc¬ 
shaped molecules via an angle and temperature depen¬ 
dent coarse-grained potential. The general purpose is to 
provide an approach that is appropriate to describe not 
only fluid phases, but also stable crystalline phases for 
these molecules. Our work builds upon a previous inves¬ 
tigation^ providing a general coarse-grained methodology 
for uniaxial discotics. In the present work, we show that 
the relevant electrostatics can be treated separately from 
the remaining interactions via a direct sum. 

In the literature there are many different approaches 
to create coarse-grained molecular models for pair inter¬ 
actions of discotics with electrostatics based on a Gay- 
Berne model^'^ or different interaction potentials.The 
simplest approach involves only one interaction site that 
inherits no orientational degree of freedom located at the 
molecules’ center-of-mass. For fluid phases at low den¬ 
sities, this kind of representation might indeed be suf¬ 
ficient, at least for molecules without complex internal 
structure.However, when considering systems at higher 
densities, the specific atomistic structure of the molecules 
yielding various degrees freedom becomes more and more 
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important. More specifically, at the other end of coarse- 
graining (before using pure quantum mechanics) are the 
atomistically-resolved models. Atomistic model studies 
for discotics have been proposed e.g. for coronene, 
benzene^^ or hexabenzororonene derivatives.^^ In our in¬ 
vestigation we make a compromise between both levels 
of detail. 

A molecular orientation vector (uniaxial) or tensor (bi¬ 
axial particle) is a frequently used example to account 
for orientational internal degrees of freedom. Such an 
approach also allows to model anisotropic shape. Gor- 
responding examples are the Gaussian-overlap^^ or the 
Gay-Berne potential^^ as well as their derivatives. 

Even more complexity can be reached by considering 
many (instead of one central) interaction sites.Fur¬ 
ther, the treatment of electrostatics can be realized 
through an electric multipole attached to the interac¬ 
tion site in a centered^’^’^ or off-centered^^arrange¬ 
ment. Moreover, by optimizing the arrangement of differ¬ 
ent electric multipoles per interaction site^^’^^ one comes 
to more and more realistic models.A quite popular 
coarse-graining strategy is given by the Martini force- 
field^^’^^ which is not based on atoms but on chemical 
building blocks. 

Here we rather aim at developing a model which has as 
few degrees of freedom as possible, while still conserving 
the uniaxiality and head-to-tail symmetry of typical disc¬ 
shaped molecules. 

Throughout the entire analysis, we focus on coronene 
molecules (see Fig. 1), for which we already developed 
an angle- and temperature dependent model for the van 
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der Waals part with desired symmetriesd Various studies 



FIG. 1. Sketch of a coronene molecule. The outer twelve 
small spheres represent hydrogen atoms while the remaining 
spheres represent carbon atoms. Equal color means equal 
distance from the center and equal charge contribution. 

suggest a broad range of applications for coronene, which 
are outlined in the following. It can be used as a building 
block for graphene nanoribbons,as a candidate for ac¬ 
tive layer compounds in photovoltaic applications^^ and 
its derivatives^^”^"^ can also be used in liquid crystal dis¬ 
plays.Understanding coronene pair interactions gives 
further insight to graphene stacking^^“^^ and growth of 
graphene.Even more, the crystal structure gives a re¬ 
lation to electronic properties, like band gaps.^^ 

We now turn to the modeling approach for coronene. 
The angle- and temperature-dependent van der Waals 
model,^ on which our investigation is based on, stems 
from a coarse-graining procedure that uses atomistic 
trajectory data from two-molecule simulations (modeled 
with the generalized Amber force-field (GAFF)^^). This 
force-field was already successfully used in growth stud¬ 
ies of another conjugated organic molecule, that is, para- 
sexyphenyl.^^ It was shown that the force field yields the 
correct phase behavior (for a more detailed discussion 
of the validation of the Amber force-field for aromatics, 
see Olivier et al.^^). Other coarse-grained potentials de¬ 
veloped in the literature^’^’^^’^^ do not take into account 
the full angle- and temperature dependence. In principle, 
the kind of coarse-graining method used in Ref. 1 could 
be applied to molecules with electrostatic contributions, 
when ab-initio simulations are used. However, ab initio 
simulations are from a computational prospective very 
time consuming compared to atomistically-detailed sim¬ 
ulations. A common compromise consists in using static 
partial charges distributed among the atoms of the un¬ 
derlying microscopic model in the coarse-graining pro¬ 
cedure.^ By implication, as shown in previous studies 
involving quantum-chemical calculations,^^ it is gener¬ 
ally not sufficient to use the static atomic partial charges 
characterizing an isolated dimer (except for selected con¬ 
figurations such as a parallel displaced one^^). How¬ 
ever, for static partial charges an angle and temperature 
dependent coarse-grained pair potential for perylene, a 
flat but non-discotic molecule, was already developed by 
Babadi et al.^ using constrained steered dynamics for 
specific configurations desired to fit an ellipsoidal soft- 


potential.^^ They also present a non-temperature depen¬ 
dent and biaxial model for coronene which implicitly in¬ 
volves static partial charges. 

An elegant approach for the electrostatics in coronene 
is presented by Obolensky et al.^^ They propose a uni¬ 
axial model for coronene consisting of concentric charged 
rings. Unfortunately, this model is not temperature de¬ 
pendent, it is based on static partial charges and the 
evaluation of the interaction potential is quite involved 
due to numerical integrations for each pair interaction, 
inconvenient for many-particle simulations. Nonetheless, 
we also consider this electrostatic contribution, which al¬ 
ready has the desired symmetries, namely uniaxiality and 
head-to-tail symmetry. 

In particular, the raw model for our investigations is 
an additive combination of the temperature-dependent 
van der Waals potential from Ref. 1 and a coarse-grained 
electrostatic potential. The entire temperature depen¬ 
dence stems from the van der Waals part alone. The 
model thus implies that the impact of slight changes in 
the charge distribution on the interatomic forces is small 
against the van der Waals forces which dominate at short 
length scales. If these temperature effects on the charge 
distribution are not important, we can focus on ground 
state charge distributions with the desired symmetry. 

We introduce two different approaches to include the 
electrostatics in a simple model. The first approach uses 
the van der Waals model and considers electrostatic in¬ 
teractions via a linear point quadrupole at the molec¬ 
ular center along the symmetry axis. This kind of ap¬ 
proach was already applied for coronene^ where an un¬ 
usual strong repulsion for planar configurations was ob¬ 
served. We address this issue in more detail later in our 
work. Nonetheless, this approach was successfully ap¬ 
plied to benzene molecules^ in the liquid phase with an 
additional dampening field for closed distances. Similar 
ideas are also used to model the interaction of clay par¬ 
ticles,which exhibit an electric double layer. 

In the second approach, we include the pure ring elec¬ 
trostatics, suggested by Obolensky et al,^"^ implicitly in 
our van der Waals model. 

We deliberately do not use any kind of point charge 
electrostatic implementations in the models since point 
charge patterns overestimate the charge localization, are 
long-ranged and do not fulfill our symmetry require¬ 
ments. 

Both of the present approaches allow a better represen¬ 
tation of the IT — IT stacking^^’^^ which is also observed in 
similar aromatic molecules, e.g. benzene dimers^^’^^ and 
hexabenzocoronene crystals. 

The remainder of this article is organized as follows. 
In section H, we introduce our models for further inves¬ 
tigation. Later on, the two model approaches are defined 
in Sec. HCl and HC2. These models are analyzed in 
Sec. HI and used for many-particle simulations in Sec. IV 
at room temperature (Sec. IV A) and beyond (Sec. IV B). 
Finally, in Sec. V we summarize our findings. 
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II. MODELS FOR PAIR INTERACTION 

A. General idea 

In the present work we assume additivity of the van der 
Waals and electrostatic contributions to the coronene- 
coronene interactions. Thus, the full pair potential is 
given through 


I/(Ra, Rb, ua, ub) = ^vdw(RA, Rb, ua, ub) 

+ ^elec.(RA, Rb, Ua, Ub), (1) 

where l/vdw aud I/giec. are the van der Waals and electro¬ 
static parts, respectively. These pair potentials are de¬ 
scribed as functions of the center of mass positions Ra, 
Rb and the axial orientations ua and ub- A configura¬ 
tion showing two discotics with these vectorial reaction 
coordinates is presented in Fig. 2. This set can be further 
reduced to four scalar reaction coordinates due to our 
symmetry requirements. These are, besides translational 
and rotational symmetry, uniaxial and head-to-tail sym¬ 
metry of the molecules (see Appendix A). The reduced 
set of reaction coordinates is important for calculating 
configuration dependent histograms (see Ref. 1). For the 
following investigations, it is more appropriate to use the 
vector-based reaction coordinates. 



FIG. 2. Sketch of a discotic dimer system with its correspond¬ 
ing reaction coordinates. 


B. Treatment of the van der Waals interaction 

A coarse-grained model describing van der Waals in¬ 
teractions has been developed in Ref. 1. Starting from an 
atomistic molecular model, involving all non-electrostatic 
interactions, we have numerically calculated the potential 
of mean force for several pair configurations and temper¬ 
atures. Then, we have parametrized the resulting curves 
using a modified Gay-Berne potential {U^qb)- The cor¬ 
responding parameters are given in Appendix D (so- 
called “vdW parametrization”). To sum up, we treat the 
van der Waals interaction using the potential 


I/vdw(R, Ua, Ub) ~ ^mGB(R, Ua, Ub) 


vdW 

parametrization 


= U:,^^s{R,ua,ub). ( 2 ) 


We stress that this parametrization is indeed 
temperature-dependent due to the rigorous coarse- 
graining procedure described in Ref. 1. Hereby 
R = Rb — Ra represents the intermolecular connecting 
vector. The mGB-potential is uniaxial in symmetry 
(suggested by the symmetry of coronene) and reads 


f^mGB(R., Ua,Ub) = 4e(R, Ua,Ub)x 



Here, i?* = (i? — (7(R, ua, ub) + rfwCTo)/(rfw o-q) repre- 
sents the reduced distance with the contact function 


cr(R,UA,UB) = 


O-Q 


1 - I (H+ + A-) + (1 - x) xt {A+ A-y 


(4) 


the well width cfq and the well depth e. The latter is 
given by 


e(R, UA, Ub) = eo [l - X^(ua • ub)^] ^ x 

X [eM(R,UA,UB)]'" • (5) 

For the coefficients in Eq. (4), we set 
A^ = ((uA ■ R) ± (ub • R))^/ (1 ± x(uA ■ Ub)) with 
the anisotropy parameters x = + 1) 

Xt = [(/^ — l)/(/D + 1)]^, where k = ctff/c^o is the quo¬ 
tient of the face-face and edge-edge contact distance. 
Regarding the well depth e [see Eq. (5)], we use the 
well-known GB formula,where the overlap factor eM 
is modified (as compared to the original definition^^) 
according to 


eM(R, Ua, Ub) — 1 “ ^ (A'^ -\r a' ) + 0 (A'^ A' )"^ 

+ ^A:(R,ua,ub). (6) 

The coefficients A'^ in Eq. (6) resemble the quantities 
A^, but incorporate the anisotropy parameter y'. The 
function 

Ar(R, Ua, Ub) = 1 - 5(ua • R) - 5(ub • R) 

- 15(ua • R)^(ub • R)^ 

-1- 2 (^(uA • Ub) - 5(ua • R) (ub • R)) (7) 

has the symmetry of a linear quadrupole-quadrupole in¬ 
teraction. 


C. Electrostatic interaction 

To get a picture of the electrostatic properties char¬ 
acterizing coronene, we start from quantum-chemical 
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results for the atomistic charge distribution given by 
Obolensky et ah (Ref. 44). As seen from Fig. 1, there are 
four different distances from the molecular center where 
the atoms are placed in a sixfold symmetry. All atoms 
with the same distance are considered to carry the same 
partial charge. We next focus on the spherical multipole 
decomposition of the electrostatic potential, that is^^ 


$(r,(/),6») = 


EE 

£=0 m=—i 


(8) 


where are the spherical harmonics depending on the 
polar angle 0 and the azimuthal angle (j). The dielectric 
constant is denoted by 5o, and the multipole moments 
are defined through 




Several symmetries in the charge distribution effectively 
lead to a reduction of the number of In particular, 

the net charge and the dipole moment of the proposed 
charge distribution are zero, i.e. (5o,o = 0 and Qi^rn = 0- 
The head-tail-symmetry {Oi = 7r/2 ^ tt — Oi) reduces 
/-values to even numbers. Summing up, our multipole 
decomposition of the electrostatic potential consists of a 
quadrupole moment {1 = 2) and higher multipoles with 

For molecules with continuous axial symmetry, we 
would only have to consider multipole moments with 
m = 0. However, the underlying atomistic model of 
coronene (see Fig. 1) does not have full (continuous) ax¬ 
ial symmetry. Nonetheless, the sixfold symmetry sug¬ 
gests that the linear quadrupole moment (Q 2 ,o) is the 
only non-zero quadrupole moment. That means the 
next non-vanishing multipole contribution after the lin¬ 
ear quadrupole is the linear hexadecapole (Q 4 ,o)- We 
next consider the quotient of the coefficients in the mul¬ 
tipole expansion corresponding to Q2,o and [see 

Eq. (8)] to evaluate the possibility of neglecting higher 
multipole terms {£ = 4 and beyond). Specifically, we 
consider 


Qe,m/{2e + l)/r^+i 

Q£,)n/(2^+ e=2,m=0,l=4,m=0 


( 10 ) 


At large distances, such as r > 2.4 nm, the absolute 
value of Q is beyond 30, i.e. the linear quadrupole ap¬ 
proach is reasonable. However, for distances smaller than 
r = 0.44 nm (a typical distance in the crystalline phase) 
the quotient Q produces values below unity. Thus, the 
quadrupole approximation is no more eligible. To some¬ 
what compensate this effect we here consider also a small 
quadrupole interaction strength to avoid an overestima¬ 
tion of the bonding. 


1. Point quadrupole approach 

In this approach the entire electrostatics is represented 
by a point quadrupole tensor. This approach is aimed 
at maintaining the electrostatics in the far field. The 
quadrupole tensor is given by^^ 

Q = ^ ® I’i - (ri)^l] > (11) 

i 

where is the charge of atom i placed at from the 
geometric center. In the eigenbasis P, the quadrupole 
tensor of a uniaxial charge distribution has the following 
symmetry 

[-1/2 0 O" 

PQP-i=Q 0 -1/2 0 . (12) 

0 0 1 

Thus, the quadrupole tensor can be written in terms of a 
linear quadrupole pointing along the mo lecula r symmetry 
axis. The prefactor in Eq. (12), Q = y/An/b Q2,o^ marks 
the quadrupole strength of the linear quadrupole. 

At this point, we can define the “point quadrupole 
model” (pq model) through the following additive for¬ 
mula 

f/pq(R, UA, ub) = UA, ub) 

where the second term on the right side describes the 
energy of one linear quadrupole in the field of another 
linear quadrupole of equal strength Q, and the function 
K is defined in Eq. (7). 

Eor simplicity, we henceforth use the dimensionless 
quadrupole moment 



The absolute is used because the sign does not influ¬ 
ence the pair potential defined in Eq. (13). The values 
of Q* resulting from density functional calculations per¬ 
formed in Ref. 44 range from 0.617 to 0.83. A well-known 
problem with a point multipole representation^’^^’^^ of 
the electrostatic interaction of extended molecules is its 
failure for closed distances, due to a spatially extended 
charge distribution and induced dipoles. In particular, 
attractive configurations become infinitely strong for van¬ 
ishing quadrupole-quadrupole distances. Indeed, we will 
later see in Sec. HI A that the latter effect leads to an 
overestimation of the binding energy for parallel dis¬ 
placed configurations. 

We here consider quadrupole strengths that are in the 
range of the previously termed reference values. Specif¬ 
ically, we use: Q* = 0,0.5,1,1.5. It is thus possible to 
gradually study the influence of the point quadrupole on 
the cohesion energy. 
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2. Charged ring approach 


The idea of this approach is to better take into account interactions for small distances. The atomic partial charges 
in the atomistic model of Obolensky et al,^^ have the same distance to the geometrical center (see Fig. 1 ) and are 
equally distributed along a ring. This makes it plausible to describe each ring as uniformly charged, thus fulfilling the 
symmetry requirements of our coarse-grained potential (uniaxiality and head-to-tail symmetry). We also note that 
considering a smeared charge distribution due to the rings is not necessarily a less accurate choice compared to an 
atomistic point charge distribution. Indeed, it seems plausible that mapping the quantum mechanical charge density 
onto atomic point charges can lead to a too strong localization. 

For simplicity, we here consider (as proposed in Ref. 44) only the two outer rings for the electrostatics. Hereby the 
hydrogen ring is charged with Qi = 12 • 0.15e and the next inner ring is oppositely charged, i.e. Q 2 = —Qi- The sum 
of all ring-ring interactions is thus given by 


rings of A rings of B 

I/ring—ring (r ^5 UA 5 ^b) — ^ ^ ^ ^ ^ring i-ring j ( 1^5 f^A? ^b) • 


Hereby a single ring-ring interaction is given by a line integral, 

Q 


Ur 


nngt-ring j(R,UA,UB) = 7 ^-^ f dr ^{Ri, sm{ 9 ), r), 
ZTTrCj J ringj 


(15) 


(16) 


where the vector r = RR + ^(ub) 


Rj cos 
Rj sin (j) 


0 


points from the center of ring i to a point on ring j. The symbol IZ 


marks a rotation matrix that rotates the x-y-plane orthogonal to ub, while Rj and (j) denote the polar coordinates of 
ring j. The electrostatic potential exerted from ring i with radius Ri on a position, described by the relative spherical 
coordinates 0 and r, is defined through 


$(i?i,sin(6»),r) = 


Qi 


2/C 


ArRi sin(0) 
r‘^+2Ri sin(0)r+it:2 


2/C 


47repm 27r 




ArRi sin(0) 

’ r‘^-2Ri sin(0)r+it:f 


-h 2rRi sm{0) -h Rf ^^2 _ 2rRi sin( 6 >) 



(17) 


where sin(^) = ||r/r x ua|| and 1C stands for the elliptic integral K. 

Direct combination of the van der Waals potential Eq. (2)] and the ring-ring potential [Eq. (15)] yields 

b^direct (r^5 f^A? ^b) ~ ^mGB(^5 f^A? ^b) T baring-ring(1^5 flA? ^b)* 


(18) 


Henceforth, we term the latter potential as “direct” 
potential. We note that t/direct is temperature- 
dependent as a consequence of the temperature- 
dependent parametrizations used for the van der Waals 
model From a computational perspective, how¬ 

ever, this kind of potential is not very suitable, at least 
not for many-particle simulations. This is, one has to 
calculate ring integrals numerically for each pair of par¬ 
ticles. Therefore, we here propose a way to effectively 
include the ring-ring interactions from above by simply 
reparametrizing the modified Gay-Berne model used for 
the van der Waals interaction [see Eq. (3)] with new pa¬ 
rameters chosen to match the “direct” potential [Eq. (18)] 
for different angular arrangements. This computationally 
more simple approach is aimed at correctly describing the 
direct potential in the near field. We name it “implicit 


electrostatics model”. Specifically, 


bbmplic. (r^5 UA? Ub) — bbnGB (r^5 UA 7 Ub) 


implic. electr. 


parametrization 

(19) 


The consequences for the long-range behavior are dis¬ 
cussed in Sec. HIB. In our implicit electrostatics 
parametrization, we fix the following parameters to the 
corresponding van der Waals values for all temperatures: 
/i = l, z^ = l, 7 = 4 and 7 ' = 4. The remaining param¬ 
eters are extracted from the direct potential [Eq. (18)] 
using the same angular configurations as those used in 
Ref. 1 for the van der Waals interaction. An overview 
of different pair configurations, including those needed to 
fit the implicit electrostatics potential [Eq. (19)] with the 
direct potential [Eq. (18)], is presented in Appendix A. 
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The anisotropy parameter y is calculated by measuring 
the face-face contact distance ctff and the edge-edge con¬ 
tact distance (Jq. The well width is calculated using ctff 
and the distance corresponding to the minimum of the 
face-face potential, i?pp^, yielding 


dw 


- (JFF 

ao (21/6 - 1)■ 


( 20 ) 


To summarize, the two parameters which determine the 
shape, are extracted from the face-face and edge-edge 
configuration. Further configurations come into play 
when we determine the remaining parameters eo, 0 
and ^ by fitting the results for Udirecti^^ ua, ub) accord¬ 
ing to Eq. (5). Specifically, our parameter fit builds 
on the four attractive wells stemming from the paral¬ 
lel weakly displaced, parallel displaced, T and edge-edge 
configuration. Finally, the specific parameters for all con¬ 
sidered temperatures are given in Appendix D (“implicit 
electrostatics parametrization”). 


III. NUMERICAL ANALYSIS OF THE MODELS 

In this section we first analyze the potential curves 
from the different electrostatic approaches described in 
Sec. lie. Second, we focus on the resulting long-range 
behavior. 


A. Potential curves 

In the following, we investigate the potential curves re¬ 
sulting from Eq. (13) (point quadrupole model), Eq. (18) 
(direct model) and Eq. (19) (implicit electrostatics 
model) for different molecular configurations at a tem¬ 
perature of T = 300 K. 

In Fig. 3, these potential curves are presented for three 
molecular pair configurations. 

Of special interest are the corresponding energy min¬ 
ima. For background information, Rapacioli et al.^^, who 
performed ground state calculations of coronene, found 
similar values for the binding energy in the face-face 
(ff) configuration (—94.9kJ/mol for tooth-to-tooth align¬ 
ment and —92.35 kJ/mol for a 30° twisted setup) and the 
parallel displaced (pd) configuration (—93.66kJ/mol). 

Similar numerical values are also found by Zhao et al.^^ 
except for the perfectly stacked face-face configuration 
(“tooth-to-tooth” setup), which is only half of the magni¬ 
tude than for twisted stacking in their calculations. Our 
present analysis shows that the direct model is character¬ 
ized by a similar binding energy for the ff-configuration 
(—50.02 kJ/mol) in comparison with the pd-configuration 
(—52.778 kJ/mol) as observed in the work of Rapacioli et 
al.^^ However, the magnitude is only about half of that 
obtained in Ref. 12. Besides the fact we consider a fi¬ 
nite temperature of 300 K and use different electrostat¬ 
ics implementations, differences also arise from the fact 



-100 ^^^^^ 
0.3 0.4 0.5 0.6 0.7 0.8 

R in nm 




1.2 1.4 1.6 1.8 

R in nm 


FIG. 3. Potential curves along the intermolecular distance 
of the direct model, point quadrupole model (with reduced 
strengths: 0, 0.5, 1) and the implicit electrostatics model at 
300 K for fixed configurations. The point quadrupole model 
with (5* = 0 coincides per definition with the van der Waals 
model. 


that our underlying molecular mechanics parameters are 
taken from the generalized AMBER force field^^ while in 
Ref. 12 parameters from van de Waal were used. A figure 
containing also ground state investigations is presented 
in Appendix B. Moreover a short analysis of the tem¬ 
perature influence on the direct model in the face-face 
configuration can be found in Appendix C. 

Anyhow, when comparing potential curves for the im¬ 
plicit electrostatics model and the direct model at 300 K, 
we find a good agreement. This is because the implicit 
electrostatics model is designed to fit the direct model in 
closed contact configurations. 

Concerning the point quadrupole model, we notice a 
strong repulsion for the ff-configuration as depicted in 
Fig. 3(a), even at the smallest considered non-vanishing 
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quadrupole strength of Q* = 0.5. This behavior was 
already observed in Ref. 7. In order to reach the poten¬ 
tial depth of the direct model in this configuration, the 
quadrupole strength should be very close to zero, which 
implies a weak far-field behavior. Moreover, even the 
smallest quadrupole strength further increases the po¬ 
tential depth of the pd-configuration [see Fig. 3(b)] away 
from its reference value. Considering with Q* = 0.5 only 
a fraction of the quadrupole strength (reference values^^ 
between Q* = 0.617 to 0.83 for various charge distri¬ 
butions; Q* = 0.963 for the proposed two-ring model), 
leads to a two times stronger potential depth for the pd- 
configuration compared to the ff-configuration of the di¬ 
rect model. This is not consistent with ground state re¬ 
sults. Nonetheless, we want to investigate in Sec. IV, 
whether a point quadrupole model is still able to predict 
the correct melting behavior of the coronene crystal. 

In Fig. 3(c), the potential curves of the edge-edge con¬ 
figuration reveal a strong influence of the electrostatics 
at all distances. This becomes clear when comparing 
^direct with the pair potential of the van der Waals model 
= ^pq(Q* = 0)- The edge-edge configuration is in 
part responsible for the binding between columnar ar¬ 
rangements in a crystal. 


B. Electrostatic long-range behavior 


The implicit electrostatics model [see Eq. (19)], which 
seems most promising so far, decays with the center-of- 
mass distance R as R~^. This means that (contrary 
to the point quadrupole model) the implicit electrostat¬ 
ics model cannot reproduce the long-range behavior ex¬ 
pected by the electrostatic interactions, which is R~^ (be¬ 
cause the quadrupole is the smallest non-vanishing multi¬ 
pole; see Sec. IIC). We note that, the ring-ring potential 
[see Eq. (15)], which is the underlying electrostatics for 
the implicit electrostatics model, reveals a correct R~^ 
decay. How does this ring-ring potential decay at short 
and intermediate distances, which are particularly rele¬ 
vant in the present context? To this end we consider the 
function 


S{R) 


d\og{Unng -ring (R, ua,ub)) 
d\og{R) 


( 21 ) 


This function is depicted in Eig. 4 for various configura¬ 
tions (see Appendix A). It turns out that S{R) is not 
approaching the far-field limit of —5 in the range of in¬ 
terest, that is, R < 2.4 nm. Therefore it is sufficient to 
model the electrostatics in our range of interest implicitly 
with a van der Waals model revealing a slightly different 
long-range behavior. A similar strategy was used to treat 
the electrostatics and the van der Waals interactions in 
DNA.^^ 




FIG. 4. Power-law exponent of the ring-ring potential for 
several molecular configurations as function of R. 


IV. MANY-PARTICLE SIMULATIONS 

In the current section we investigate a many-particle 
system of coronene molecules using molecular dynamics 
at constant temperature (T) and pressure (P). First, we 
investigate the equilibrium structural properties of the 
unit cell at room temperature. Second, we analyze how 
certain structural parameters behave upon heating up 
the crystal until it melts. Each simulation run starts from 
lattice-like initial conditions (described below) and leads 
to a relaxation of particle positions and box-vectors. To 
affect box lengths as well as box-angles we use anisotropic 
pressure coupling. In particular, the temperature and 
pressure coupling was realized with the Berendsen weak 
coupling scheme,also used in our previous work.^ Ro¬ 
tational dynamics was solved in a similar fashion, as pro¬ 
posed by Eincham et al.^^ with separate temperature con¬ 
trol. 

The coronene molecules are modeled with either the 
point quadrupole model (see Sec. IICl) or the implicit 
electrostatics model (see Sec. IIC 2). Our many-particle 
system for all simulations consists of 576 molecules and is 
initialized with the following arrangement of unit cells in 
terms of box-vectors: Li = 4 a, L2 = 12 b and L3 = 6 c. 
Hereby a, b and c denote the right-handed unit cell vec¬ 
tors, which together with the molecular arrangement and 
orientation in the unit cell are taken from Ref. 60. The 
molecules are aligned in a so-called herringbone pattern, 
which resembles tractor traces when seen from the side 
(see Eig. 5). 

The pressure is set to 1 bar and the compressibility to 
2.25 • 10“^bar~^. Relaxation time constants involved in 
T- and P-control are set to 2 ps. The equilibration times 
range between 3.5 and 100 ns at a time step of 10 fs. A 
force-shifted cutoff was used for the pair forces and set 
to 2.4 nm. To test this approach, we calculated the total 
electrostatic energy both, with the cutoff and by using 
the full Ewald sum for quadrupoles (for the explicit ex¬ 
pression, see Ref. 61). We found that the electrostatic 
energy from the cutoff procedure differs from the corre- 

















(a) 



(b) 


FIG. 5. (a) Unit cell with molecules A and B (the red lines 
represent the side view of the discotic molecules), (b) Snap¬ 
shot of a real crystalline configuration. 


spending energy using the Ewald summation technique 
by only half a per mille. This justifies the simpler cut-off 
procedure. 


A. Bulk crystal at room temperature 


In the following, we investigate the equilibrium struc¬ 
ture of the bulk crystal for the point quadrupole model 
(see Sec. IIC 1 ; quadrupole strengths: Q* = 0 , 0.5, 1 ,1.5) 
and the implicit electrostatics model (see Sec. IIC 2 ) us¬ 
ing trajectory data from the Molecular Dynamics simu¬ 
lations. One important criteria to judge the performance 
of the different models are the unit cell parameters listed 
in Tab. I. 


TABLE I. Variation of unit cell parameters for the different pair potentials at room temperature. The box lengths of the cell 
are denoted with a, b and c defined through the length of the unit cell vectors a, b and c. The angles between the unit cell 
vectors are denoted with a = k(b, c), /3 = Z(a, c) and 7 = Z(a, b). The volume and the cohesion energy per particle are listed 
as well. 


Model 

a(A) 

6(A) 

c(A) 

a(°) 

/3(°) 

7(°) 

y(A®) 

-E'coh(eV) 

pq model (Q*=0) / vdW model 

19.6444 

3.6282 

9.1430 

90.02 

117.74 

89.98 

576.78 

1.5533 

pq model (Q*=0.5) 

17.0045 

4.2382 

9.2411 

90.00 

113.66 

90.00 

610.04 

1.7153 

pq model {Q* = 1) 

17.1237 

4.0178 

9.7228 

90.00 

116.52 

90.00 

598.53 

2.9882 

pq model (Q*=1.5) 

17.7104 

3.7246 

9.9783 

89.62 

122.77 

86.83 

552.08 

6.4402 

impl. electr. model 

17.6142 

4.5524 

9.5722 

90.00 

112.91 

90.01 

707.03 

0.9457 

experiment®^ 

16.094(9) 

4.690(3) 

10.049(8) 

90 

110.79(2) 

90 

709.9(8) 

/ 

experiment®®’®® 

16.119 

4.702 

10.102 

90 

110.9 

90 

717.1 

1.39 - 1.54 

database®^ (-^coh from DFT-calculations®®) 

16.11 

4.70 

10.10 

90 

110.9 

90 

714.43 

1.378-1.783 


Inspecting the values for a, 5, c, a, 7 (for descrip¬ 
tion see caption of Tab. I), we find very good agree¬ 
ment of the implicit electrostatics model and the point 
quadrupole model for weak strength (Q* = 0.5) with the 
corresponding experimental data. The remaining mod¬ 
els (Q* 7 ^ 0.5) reveal slight deviations in a, 6 , c, and p. 
Nonetheless, the crystal structures are in all cases mono¬ 
clinic (a = 7 = 90°). Concerning the cohesion energy pre¬ 
dicted by the point quadrupole models, we note a signif¬ 
icant increase with the quadrupole strength. This stems 


from the overestimation of the binding energy for par¬ 
allel displaced configurations (see Eig. 3). The cohesion 
energy of the implicit electrostatics model is somewhat 
underestimated. 

We further want to analyze the orientation of the cor¬ 
responding molecules (named A and B) in the unit cell. 
Therefor, we introduce the herringbone angle 0, which 
marks the angle between both molecular orientation axes 
that point “upwards” along b (see Eig. 5), and is defined 
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as follows 

(j) = arccos [(sgn(uA * b)uA) • (sgn(uB • b)uB)] • (22) 

In Table II all molecular angles with respect to the 
body fixed unit cell are displayed. Specifically, we present 
the enclosed angles of the molecule A’s axis with each 
unit cell axis as performed in Ref. 60. 


TABLE II. The angles of molecule A with the a, b and its 
perpendicular axis a x b are denoted by . The 

notation is taken from Ref. 60. 


Model 


x%n 


0 

3 

0(°) 

pq model {Q* 

= 0) / vdW model 

95.65 

12.20 

79.24 

24.38 

pq model {Q* 

= 0.5) 

129.69 

39.86 

93.12 

79.72 

pq model {Q* 

= 1) 

129.83 

40.66 

96.85 

81.32 

pq model {Q* 

= 1.5) 

121.90 

40.79 

108.47 

81.54 

impl. electr. model 

130.23 

40.31 

92.30 

80.63 

ideal crystal®^ 


133.7 

43.7 

89.6 

87.35 


As expected, we observe a nematic phase for the van 
der Waals model (point quadrupole model with a zero 
quadrupole strength) reflected by a small value for the 
enclosed angle of the molecule A’s axis with the b-axis, 
called and a small herringbone angle (j). This finding 
is in accord with our previous work.^ The alignment of 
the molecules with respect to all unit cell axes {Xn^ 
uj^) is for both, the implicit electrostatics model and the 
point quadrupole model with weak strength (Q* = 0.5), 
in good agreement with the experimental values. For 
all investigations (except for the nematic phase of the 
model without electrostatics), no significant change for 
the distance of neighboring columns was observed (cor¬ 
responding distances are 8.7 A and 9.8 A). 


B. Melting of bulk crystal 

We now turn to the investigation of structural changes 
of the bulk crystal [see Fig. 5(b)] upon heating up the 
system. Hereby, we expect the crystal to melt. The tem¬ 
perature range considered covers 300 K < T < 1500 K. 
After every 100 K a different isothermal-isobaric simula¬ 
tion run is performed. 

We further investigate 0 as a function of temperature 
for the different models as presented in Fig. 6 (a). The 
melting from crystalline to isotropic phases is also re¬ 
flected in Fig. 6 (b) showing the third root of the volume 
as a function of the temperature. 

When considering the reference temperatures of bulk 
coronene with respect to melting (710.5 K)^^ and boiling 
(798K)^^, we recognize that coronene is liquid only in a 
narrow temperature band. In the implicit electrostatics 
model we observe a decay of the herringbone angle be¬ 
tween 600 and 700 K towards zero indicating isotropic ori¬ 
entation. At 700 K the volume is significantly increased 



pq model (Q* = 0 ) — 
pq model ((5* = 0.5) — 
pq model (Q* = 1 ) — 
pq model {Q* = 1.5) -r- 
impl. electr. model — 


FIG. 6. (a) Herringbone angle between neighboring discs and 
(b) volume as function of the temperature. 


in comparison to 600 K, but highly raises with increas¬ 
ing temperature. We take these as indications of a liq¬ 
uid phase appearing at around 700 K and a gas phase 
at temperatures bigger than 800 K. On the contrary to 
the implicit electrostatics model, the crystal structure for 
the point quadrupole model melts at significantly higher 
temperatures for all quadrupole strengths. It seems ob¬ 
vious that the melting temperature increases with the 
cohesion energy, as discussed before (see Sec. IV A). 

Finally, we want to analyze the crystalline order paral¬ 
lel to the plane spanned by Li and L 3 with the following 
correlation function^^ 

9\\{h) = /(/i,n-Rjfe)\ , (23) 

\ II ^ j fe>j / 

where p is the density and f{x^y) equals unity for 
y ^ (x — y,x+y], otherwise / = 0 (with A being the 
bin size). The volume appearing in Eq. (23) is defined as 
I^l = A |Li X Lsj. For the implicit electrostatics model, 
we observe a crystalline order [see Fig. 7(a)], that con¬ 
tinuously vanishes with rising temperature. In contrast, 
the point quadrupole models do not exhibit this behavior 
[see Fig. 7(b)-(e)]. At a quadrupole strength of Q* = 1.5 
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we do not see a positional order any more, while the ori¬ 
entational order still exists (see Fig. 6). 

To sum up, the transition temperature of the implicit 
electrostatics model can be suitably reproduced. Clearly 
it is not in full accordance with the experimental refer¬ 
ence values. Nonetheless, this approach marks a way to 
treat the relevant electrostatics in a rather simple model 
eligible for large-scale computer simulations. 

V. CONCLUSION 

In this study we presented different approaches to 
model the full pair interaction (with electrostatics) be¬ 
tween coronene molecules, starting from a previously in¬ 
troduced van der Waals-like of model (see Sec. IIB) that 
neglected electrostatics. Hereby, coronene serves as an 
exemplary discotic molecule with head-to-tail symmetry. 
In the first approach (see Sec. IICl), we aimed at ex¬ 
tending this van der Waals model by a point quadrupole 
attached at the particles’ center of mass along the par¬ 
ticles’ symmetry axis. In the second approach (see 
Sec. IIC 2) the van der Waals model was extended via 
electrostatic interactions stemming from charged rings. 
These ring interactions were used to reparametrize the 
van der Waals model into the so-called implicit elec¬ 
trostatics model. In both electrostatic approaches, we 
treated the relevant electrostatics separate from the van 
der Waals interactions via a direct sum. 

To assess the quality of modeling, unit cell parame¬ 
ters and structural quantities of the coronene bulk crys¬ 
tal were calculated with constant pressure and tempera¬ 
ture Molecular Dynamics. These results have been com¬ 
pared to their experimentally determined counterparts 
(see Sec. IV). 

Based on our data, we can conclude that the simple 
point quadrupole approach, although yielding the correct 
long-range decay gives unreliable results, even when the 
quadrupole strength is reduced to a value related to a 
reasonable cohesion energy. 

Moreover, for all considered quadrupole strengths, we 
observed melting temperatures far beyond the experi¬ 
mental values. Nevertheless, a small point quadrupole 


leads to a stabilization of the herringbone structure due 
to an energetic preference. In contrast, our second ap¬ 
proach involving the implicit electrostatics model is able 
to stabilize the herringbone structure up to a melting 
point similar to experimental values. With this second 
model, we also encountered a liquid phase as observed 
in the experiments. As a conclusion, the implicit elec¬ 
trostatics model seems superior in reproducing structure 
and melting. It is also more convenient from a com¬ 
putational perspective since it is just a reparametrized 
Gay-Berne-like model. The crucial point for the success 
of the latter approach consists of choosing an extended 
charge distribution rather than a single point multipole. 
Hereby the characteristics of the electrostatic potential 
in the near field are reproduced in a great extent. Still, 
it would be interesting to investigate the influence of an 
atomic point charge distribution on the structural and 
melting properties of the coronene crystal. 

Besides the temperature-independent electrostatics 
and the ring-charge assumption, a further underlying as¬ 
sumption in our work is the pairwise additivity of the 
many-molecule interactions. The implications of this as¬ 
sumption were already investigated for crystalline ben¬ 
zene,which is similar to coronene. Further, it is worth 
mentioning that our temperature-dependent approach 
does not take into account formation and breaking of co¬ 
valent bonds, which may be an important processes for 
very high temperatures and pressures. 

For future work, it would be rather interesting to use 
ab-initio simulations to calculate effective pair potentials 
for coronene and compare them with our corresponding 
results for different temperatures. Finally, we propose 
applicability of our presented method to a similar class of 
discotic molecules, e.g. benzene or hexabenzocoronene. 
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Appendix A: Dimer configurations 


In Table HI we present important examples of dimer configurations in terms of four reaction coordinates. The 
reaction coordinates R, a,5 and c are related to the molecular center of mass positions Ra, Rb and the orientations 
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' ■ ■ ' I ' ' ' ' I ■ ' ' ' I 

(a) implicit electr. 
model 


■ ' ' ' I ' ' ' ' I ' ' ' ■ I ■ ' 

(b) pq model ((5*=0) 


' ■ ■ ■ I ' ' ' ' I ' ' ' ' I ' ■ ■ ' 

(c) pq model (Q* = 0.5) 


^2 



0 0.2 0.4 0.6 0 0.2 0.4 0.6 0 0.2 0.4 0.6 

h in nm h in nm h in nm 



300A: - 

ioooa: - 

400A: - 

iiooa: - 

500A: - 

1200A: - 

600A: - 

1300a: -- 

700A: - 

1400A: - 

800A: - 

900A: - 

1500A: - 


FIG. 7. Temperature influence of pair correlation function along the Li-La-plane for the implicit electrostatics model and the 
point quadrupole model for various strengths. 


ua and up as follows 


R = 11 Ftp ~ IF. 


a = 


b = 


Ua R 


Up R 


■All , 


with R = 


Rp — R^ 

||Rb - 


c = sgn(uA • R) sgn(up • R) ua • up. 


(Al) 

(A2) 

(A3) 

(A4) 


These coordinates describe a molecular pair in the body-fixed frame, providing that the particles are uniaxial and 
have head-to-tail symmetry. In addition a, b and c equally treat chiral dimer configurations. 


TABLE III. This table summarizes the interesting configurations with their corresponding values of the reaction coordinates 
according to Ref. 1. 



face-face 

H 

a b c 

parallel weakly 

displaced 
a b c 

parallel displaced 

// 

a b c 

T 

•• 

a b c 

herringbone 

a be 

V 

1/ 

a b c 

edge-edge 

a b c 

cross 

o- 

a b c 

value 

11 1 

0.94 0.94 1 

I/V 2 I/V 2 1 

0 10 

0.8767 0.481 0 

1 I/V 2 I/V 2 

0 0 1 

0 0 0 
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Appendix B: Ground state curves 


In order to compare our results with ground state (GS) results in the literature, we consider in this Appendix a 
further model, which combines the van der Waals potential pertaining to the atomistic level (using generalized 

AMBER force field"^^ without partial charges) at T = OK and our ring-ring potential from Eq. (15). Specifically, 




direct 




vdW 


U ring—ring- 


(Bl) 


We term this potential “ground state direct model”. In Eig. 8 we present the ground state potential curves 
(without electrostatics) and [with ring-ring electrostatics from Eq. (15)] using various configurations suitable 

for comparison with literature results.In addition, we add the potential curves showing the ground state analogon 
of the point quadrupole model, that is 

qT = fiA, "b) + ^ ^ ■ K{K, ua, ub). (B2) 



0.7 0.8 0.9 1 1.1 

R in nm 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


i?ll in nm 



R in nm 


FIG. 8. Potential curves of the ground state models ^^d Ihe direct model Edirect, the implicit electrostatics model 

Eimpiic. and the point quadrupole model (reduced strengths: 0.5, 1, 1.5) at T = 300 K. (a) face-face configuration (with 
tooth-to-tooth setup); (b) parallel displaced configurations for a layer distance of 0.332nm; (c) T configuration; (d) edge-edge 
configuration. 


We recognize that characterized by similar 

values of the binding energy for the ff-configuration in 
Eig. 8(a) (—58.31 kJ/mol with tooth-to-tooth setup) and 
the pd-configuration in Eig. 8(b) (—61.37 kJ/mol), con¬ 
sistent to what was observed in the work of Rapacioli 
et al^^ (ff-configuration: —94.9kJ/mol, pd-configuration: 


—93.66 kJ/mol). However, the magnitudes of 
only about two thirds of that obtained in Ref. 12. These 
differences stem from different charge distributions used 
by us and in Ref. 12, and from the fact that our atomistic 
parameters are taken from the generalized AMBER force 
field^^ while in Ref. 12 parameters from van de Waal were 
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used. Still, we conclude that the curves of i*^veal 

far more lower binding energies than the 300 K curves 
of t/direct ^^e thus closcr to the mentioned literature 
values. 

For the sake of completeness, we also present curves for 
the t- and edge-edge configuration in Fig. 8(c) and 8(d). 
The general impact of the electrostatics is reflected by 
the difference between quite pro¬ 

nounced in all configurations. As mentioned in Sec. Ill A 
the point quadrupole curves strongly affect the binding 
energies for all configurations. 


We can state that the higher the temperature, the flat¬ 
ter the potential minimum and the larger the binding dis¬ 
tance. The temperature dependence mainly stems from 
molecular bending modes as discussed in Ref. 1. Its de¬ 
pendence due to axial averaging is rather small, reflected 
by the small difference of both ground state potentials. 
It is also worthwhile to mention that there is no influence 
of temperature at around R = 0.42 nm, because at this 
point all potential curves intersect. 


Appendix C: Temperature dependence 


The temperature dependence of the face-face configu¬ 
ration is shown in Fig. 9 for the direct model and the 
ground state direct model [see Eq. (BI)] for two different 
configurations (tooth-to-tooth and 30°-twisted). 



FIG. 9. Potential curves for the face-face configuration using 
the direct model at 300 K, 800 K and 1500 K. In addition two 
ground state configurations are added. 


Appendix D: Parametrizations 


The following tables summarize the parameters used in Eq. (3) for the van der Waals model (Tab. IV) and the 
implicit electrostatics model (Tab. V) at different temperatures T. 
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TABLE IV. Van der Waals parametrization 


T(K) 

K 

x' 

g 

V 

cro(nm) 

Co (kJ / mol) 

d^ 

7 

1 

e 


300 

0.2885 

-0.7895 

1 

1 

1.0529 

6.5481 

0.3884 

4 

4 

0.1592 

-0.1967 

400 

0.2892 

-0.7916 

1 

1 

1.0603 

6.3553 

0.3869 

4 

4 

0.2019 

-0.1984 

500 

0.2899 

-0.7939 

1 

1 

1.0678 

6.1614 

0.3854 

4 

4 

0.2478 

-0.2002 

600 

0.2905 

-0.7962 

1 

1 

1.0752 

5.9666 

0.3839 

4 

4 

0.2971 

-0.2022 

700 

0.2912 

-0.7986 

1 

1 

1.0826 

5.7707 

0.3825 

4 

4 

0.3503 

-0.2043 

800 

0.2919 

-0.8012 

1 

1 

1.0900 

5.5739 

0.3811 

4 

4 

0.4078 

-0.2065 

900 

0.2916 

-0.8024 

1 

1 

1.0906 

5.4425 

0.3901 

4 

4 

0.4426 

-0.2063 

1000 

0.2914 

-0.8037 

1 

1 

1.0912 

5.3113 

0.3992 

4 

4 

0.4790 

-0.2060 

1100 

0.2911 

-0.8051 

1 

1 

1.0918 

5.1802 

0.4082 

4 

4 

0.5171 

-0.2057 

1200 

0.2909 

-0.8064 

1 

1 

1.0925 

5.0494 

0.4172 

4 

4 

0.5571 

-0.2055 

1300 

0.2907 

-0.8079 

1 

1 

1.0931 

4.9187 

0.4262 

4 

4 

0.5991 

-0.2052 

1400 

0.2904 

-0.8094 

1 

1 

1.0937 

4.7882 

0.4352 

4 

4 

0.6434 

-0.2049 

1500 

0.2902 

-0.8109 

1 

1 

1.0943 

4.6579 

0.4442 

4 

4 

0.6900 

-0.2046 


TABLE V. Implicit electrostatics parametrization 


T(K) 

K, 

x' 

g 

V 

ao(nm) 

eo(kJ/mol) 

dw 

7 

l' 

9 


300 

0.2856 

-0.7639 

1 

1 

1.1026 

5.3125 

0.3440 

4 

4 

1.4825 

-0.3123 

400 

0.2867 

-0.7635 

1 

1 

1.1088 

5.2244 

0.3406 

4 

4 

1.5003 

0.3120 

500 

0.2877 

-0.7631 

1 

1 

1.1149 

5.1354 

0.3373 

4 

4 

1.5196 

-0.3117 

600 

0.2888 

-0.7627 

1 

1 

1.1211 

5.0456 

0.3340 

4 

4 

1.5403 

-0.3114 

700 

0.2898 

-0.7623 

1 

1 

1.1273 

4.9548 

0.3307 

4 

4 

1.5626 

-0.3111 

800 

0.2909 

-0.7618 

1 

1 

1.1335 

4.8633 

0.3275 

4 

4 

1.5865 

-0.3107 

900 

0.2905 

-0.7622 

1 

1 

1.1362 

4.7541 

0.3378 

4 

4 

1.6406 

-0.3111 

1000 

0.2902 

-0.7626 

1 

1 

1.1390 

4.6452 

0.3480 

4 

4 

1.6969 

-0.3115 

1100 

0.2899 

-0.7631 

1 

1 

1.1418 

4.5365 

0.3582 

4 

4 

1.7558 

-0.3119 

1200 

0.2896 

-0.7635 

1 

1 

1.1446 

4.4280 

0.3683 

4 

4 

1.8173 

-0.3123 

1300 

0.2892 

-0.7639 

1 

1 

1.1474 

4.3197 

0.3784 

4 

4 

1.8816 

-0.3127 

1400 

0.2889 

-0.7644 

1 

1 

1.1501 

4.2117 

0.3885 

4 

4 

1.9491 

-0.3132 

1500 

0.2886 

-0.7649 

1 

1 

1.1529 

4.1038 

0.3984 

4 

4 

2.0198 

-0.3137 
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